home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / CUGUK / C005.ZIP / MATH.C < prev    next >
Text File  |  1990-01-19  |  13KB  |  533 lines
  1. /********************************************************************
  2.  * C Users Group (U.K) C Source Code Library File CUGLIB.005        *
  3.  * Inquiries to: M. Houston, 36 Whetstone Clo. Farquhar Rd.         *
  4.  * Edgbaston, Birmingham B15 2QN ENGLAND                *
  5.  ********************************************************************
  6.  * File name: math.c
  7.  * Program name: library modules only
  8.  * Source of file: The Public Domain Software Library.
  9.  * Purpose: A collection of common maths functions.
  10.  * Changes: <who what when & why major changes have been made>      
  11.  ********************************************************************/
  12.  
  13. /*****************************************************************************
  14.  *                                         *
  15.  * MATH.C  Mathematical subroutines for use with Lattice or Microsoft-         *
  16.  *       compiler.                                 *
  17.  *       By Max R. Darsteler                             *
  18.  *                                         *
  19.  *****************************************************************************
  20.  */
  21.  
  22. /*****************************************************************************
  23.  *  This file contains the combined source files of the most
  24.  *  important mathe- matical formulas written for the Lattice C Compiler.
  25.  *  It makes no use of the 8087.  Some of the functions make use of the
  26.  *  internal binary representation of double precision numbers.  As the
  27.  *  internal data representation may vary between different compilers,
  28.  *  watch out for errors e.g.  in the power function.  The function values
  29.  *  are calculated through polynominal or rational approximations.  The
  30.  *  results are for allmost all functions precise for the first 9 to 10
  31.  *  decimal digits.  For still higher precision or checking of values use
  32.  *  functions from file mathf.c (not available here), which contains
  33.  *  Taylor's expansions series.  As this requires many iterations with
  34.  *  multiplications and divisions, the calculation time is too long for use
  35.  *  in regular programs.  In order to use these functions, a program has to
  36.  *  include the header file "math.h".  Furthermore, either the compiled
  37.  *  version (math.obj) or a library with the compiled functions has to be
  38.  *  linked to the program.  As exemples see included testprograms
  39.  *  "tstsin.c" and "tstpow.c" To compile:  "mcb tstsin <CR>", to link:"link
  40.  *  c+%1+math.obj, %1, ,mc <CR>".  I recommend to compile the source of
  41.  *  each function seperately with your personal compiler and use the MS-DOS
  42.  *  librarian to build a library of just the functions you need.  Some of
  43.  *  the more esotheric functions used on "UNIX" are such as Bessel
  44.  *  functions or the gamma function are not implemented.  Evidently, this
  45.  *  is not a professional implementation.
  46.  *
  47.  *   Rockville, 11/27/83      Max R. Drsteler
  48.  *                   12405 Village Sq. Terrace
  49.  *                   Rockville Md. 20852    H 301 984-1168
  50.  ******************************************************************************/
  51.  
  52. /* SQRT.C
  53.  * Calculates square root
  54.  * Precision: 11.05 E-10
  55.  * Range: 0 to e308
  56.  * Author: Max R. Drsteler
  57.  */
  58. double sqrt(x)
  59. double x;
  60. {
  61.     typedef union {           /* Makes use of internal data  */
  62.         double d;          /* representation */
  63.         unsigned u[4];
  64.     } DBL;
  65.     double y, re, p, q;
  66.     unsigned ex;        /* exponens*/
  67.     DBL *xp;
  68.  
  69.     xp = (DBL *)&x;
  70.     ex = xp->u[3] & ~0100017; /* save exponens */
  71.     re = ex & 020 ? 1.4142135623730950488 : 1.0;
  72.     ex = ex - 037740 >> 1;
  73.     ex &= ~0100017;
  74.     xp->u[3] &= ~0177760;      /* erase exponens and sign*/
  75.     xp->u[3] |= 037740;      /* arrange for mantissa in range 0.5 to 1.0 */
  76.     if (xp->d < .7071067812) {  /* multiply by sqrt(2) if mantissa < 0.7 */
  77.         xp->d *= 1.4142135623730950488;
  78.         if (re > 1.0) re = 1.189207115002721;
  79.         else re = 1.0 / 1.18920711500271;
  80.     }
  81.     p =           .54525387389085 +   /* Polynomial approximation */
  82.         xp->d * (13.65944682358639 +   /* from: COMPUTER APPROXIMATIONS*/
  83.         xp->d * (27.090122712655   +   /* Hart, J.F. et al. 1968 */
  84.         xp->d *   6.43107724139784));
  85.     q =          4.17003771413707 +
  86.         xp->d * (24.84175210765124 +
  87.         xp->d * (17.71411083016702 +
  88.         xp->d ));
  89.     xp->d = p / q;
  90.     xp->u[3] = xp->u[3] + ex & ~0100000;
  91.     xp->d *= re;
  92.     return(xp->d);
  93. }
  94.  
  95.  
  96. /* EXP2.C
  97.  * Calculates exponens to base 2 using polynomial approximations
  98.  * Precision: 10E-10
  99.  * Range:
  100.  * Author: Max R. Drsteler
  101.  */
  102. double exp2(x)
  103. double x;
  104. {
  105.     typedef union {
  106.         double d;
  107.         unsigned u[4];
  108.     } DBL;
  109.     double y, x2, p, q, re;
  110.     int ix;
  111.     DBL *xp, *yp;
  112.  
  113.     xp = (DBL *)&x;
  114.     y = 0.0;
  115.     yp = (DBL *)&y;
  116.     if(xp->d > 1023.0 || xp->d < -1023.0) return (1E307);
  117.     ix = (int) xp->d;
  118.     yp->u[3] += ix + 1023 << 4;
  119.     yp->u[3] &= ~0100017;
  120.     if ((xp->d -= (double) ix) == 0.0) return(yp->d);
  121.     if (xp->d < 0.0) {
  122.           yp->u[3] -= 1 << 4;
  123.           yp->u[3] &= ~0100017;
  124.           xp->d++;
  125.     }
  126.     if (xp->d >= 0.5) {  /* adjust to range 0-0.5 */
  127.         xp->d -= 0.5;
  128.         re = 1.41421356237309504880;
  129.     }
  130.     else re = 1.0;
  131.     x2 = xp->d * xp->d;
  132.     p = xp->d * (7.2152891511493 +
  133.            x2 *  0.0576900723731);
  134.     q =        20.8189237930062 +
  135.            x2 ;
  136.     xp->d = (q + p) / (q - p);
  137.     yp->d *= re * xp->d;
  138.     return(yp->d);
  139. }
  140.  
  141.  
  142. /* POW.C
  143.  * Calculates y^x
  144.  * Precision:
  145.  * Range: o to big for y, +/-1023 for x
  146.  * Author: Max R. Drsteler, 10/2/83
  147.  */
  148.  
  149. double pow(y,x)
  150. double y, x;
  151. {
  152.     typedef union {
  153.         double d;
  154.         unsigned u[4];
  155.     } DBL;
  156.     double z, w, p, p2, q, re;
  157.     unsigned ex;        /* exponens*/
  158.     int iz;
  159.     DBL *yp, *zp, *wp;
  160.  
  161.     yp = (DBL *)&y;
  162.     if (yp->d <= 0.0) y = -y;
  163.     z = 0.0;
  164.     zp = (DBL *)&z;
  165.     zp->u[3] = yp->u[3] & ~0100017; /* save exponens */
  166.     iz = (zp->u[3] >> 4)-1023;
  167.     if ((yp->d - zp->d) == 0.0)
  168.         z = (double)iz;
  169.     else {
  170.         yp->u[3] -= ++iz << 4; /* arrange for range 0.5 to 0.99999999999 */
  171.         yp->d *= 1.4142135623730950488;  /* shift for 1/sqrt(2) to sqrt(2) */
  172.         p = (yp->d - 1.0) / (yp->d + 1.0);
  173.         p2 = p * p;
  174.         z = p  * (2.000000000046727 +    /* Polynomial approximation */
  175.             p2 * (0.666666635059382 +    /* from: COMPUTER APPROXIMATIONS*/
  176.             p2 * (0.4000059794795   +    /* Hart, J.F. et al. 1968 */
  177.             p2 * (0.28525381498     +
  178.             p2 *  0.2376245609 ))));
  179.         z = z * 1.442695040888634 + (double)iz    - 0.5;
  180.     }
  181.     z *= x;
  182.     w = 0.0;
  183.     wp = (DBL *)&w;
  184.     if(zp->d > 1023.0 || zp->d < -1023.0) return (1E307);
  185.     iz = (int) zp->d;
  186.     wp->u[3] += iz + 1023 << 4;
  187.     wp->u[3] &= ~0100017;
  188.     if ((zp->d -= (double) iz) == 0.0) return(wp->d);
  189.     while (zp->d < 0.0) {
  190.           wp->u[3] -= 1 << 4;
  191.           wp->u[3] &= ~0100017;
  192.           zp->d++;
  193.     }
  194.     if (zp->d >= 0.5) {  /* adjust to range 0-0.5 */
  195.         zp->d -= 0.5;
  196.         re = 1.41421356237309504880;
  197.     }
  198.     else re = 1.0;
  199.     p2 = zp->d * zp->d;
  200.     p = zp->d * (7.2152891511493 +
  201.            p2 *  0.0576900723731);
  202.     q =        20.8189237930062 +
  203.            p2 ;
  204.     zp->d = re * wp->d * (q + p) / (q - p);
  205.     return(zp->d);
  206. }
  207.  
  208.  
  209. /* FMOD.C
  210.  * Returns the number f such that x = i*y + f. i is an integer, and
  211.  * 0 <= f < y.
  212.  */
  213. double fmod(x, y)
  214. double x, y;
  215. {
  216.     double zint, z;
  217.     int i;
  218.  
  219.     z = x / y;
  220.     zint = 0.0;
  221.     while (z > 32768.0) {
  222.         zint += 32768.0;
  223.         z -= 32768;
  224.     }
  225.     while (z < -32768.0) {
  226.         zint -= 32768.0;
  227.         z += 32768.0;
  228.     }
  229.     i = (int) z;
  230.     zint += (double) i;
  231.     return( x - zint * y);
  232. }
  233.  
  234. /*ASIN.C
  235.  *Calculates arcsin(x)
  236.  *Range: 0 <= x <= 1
  237.  *Precision: +/- .000,000,02
  238.  *Header: math.h
  239.  *Author: Max R. Drsteler
  240.  */
  241. extern double sqrt();
  242.  
  243. double asin(x)
  244. double x;
  245.  
  246. {
  247.     double y;
  248.     int sign;
  249.  
  250.     if (x > 1.0 || x < -1.0) exit(1);
  251.     sign = 0;
  252.     if (x < 0) {
  253.          sign = 1;
  254.          x = -x;
  255.     }
  256.     y = ((((((-.0012624911    * x
  257.          + .0066700901) * x
  258.          - .0170881256) * x
  259.          + .0308918810) * x
  260.          - .0501743046) * x
  261.          + .0889789874) * x
  262.          - .2145988016) * x
  263.          +1.5707963050;
  264.     y = 1.57079632679 - sqrt(1.0 - x) * y;
  265.     if (sign) y = -y;
  266.     return(y);
  267. }
  268.  
  269. /* SIN.C
  270.  * Calculates sin(x), angle x must be in rad.
  271.  * Range: -pi/2 <= x <= pi/2
  272.  * Precision: +/- .000,000,005
  273.  * Header: math.h
  274.  * Author: Max R. Drsteler
  275.  */
  276.  
  277. double sin(x)
  278. double x;
  279. {
  280.     double xi, y, q, q2;
  281.     int sign;
  282.  
  283.     xi = x; sign = 1;
  284.     while (xi < -1.57079632679489661923) xi += 6.28318530717958647692;
  285.     while (xi > 4.71238898038468985769) xi -= 6.28318530717958647692;
  286.     if (xi > 1.57079632679489661923) {
  287.         xi -= 3.141592653589793238462643;
  288.         sign = -1;
  289.     }
  290.     q = xi / 1.57079632679; q2 = q * q;
  291.     y = ((((.00015148419  * q2
  292.           - .00467376557) * q2
  293.           + .07968967928) * q2
  294.           - .64596371106) * q2
  295.           +1.57079631847) * q;
  296.     return(sign < 0? -y : y);
  297. }
  298.  
  299.  
  300. /* LOG.C
  301.  * Calculates natural logarithmus
  302.  * Precision: 11.56 E-10
  303.  * Range: 0 to e308
  304.  * Author: Max R. Drsteler
  305.  */
  306. double log(x)
  307. double x;
  308. {
  309.     typedef union {
  310.         double d;
  311.         unsigned u[4];
  312.     } DBL;
  313.     double y, z, z2, p;
  314.     unsigned ex;        /* exponens*/
  315.     int ix;
  316.     DBL *xp, *yp;
  317.  
  318.     xp = (DBL *)&x;
  319.     if (xp->d <= 0.0) return(y);
  320.     y = 0.0;
  321.     yp = (DBL *)&y;
  322.     yp->u[3] = xp->u[3] & ~0100017; /* save exponens */
  323.     ix = (yp->u[3] >> 4)-1023;
  324.     if ((xp->d - yp->d) == 0.0) return( .693147180559945 * (double)ix);
  325.     xp->u[3] -= ++ix << 4; /* arrange for range 0.5 to 0.99999999999 */
  326.     xp->d *= 1.4142135623730950488;  /* shift for 1/sqrt(2) to sqrt(2) */
  327.     z = (xp->d - 1.0) / (xp->d + 1.0);
  328.     z2 = z * z;
  329.     y = z  * (2.000000000046727 +    /* Polynomial approximation */
  330.         z2 * (0.666666635059382 +    /* from: COMPUTER APPROXIMATIONS*/
  331.         z2 * (0.4000059794795   +    /* Hart, J.F. et al. 1968 */
  332.         z2 * (0.28525381498     +
  333.         z2 *  0.2376245609 ))));
  334.     y = y + .693147180559945 * ((double)ix    - 0.5);
  335.     return(yp->d);
  336. }
  337.  
  338. /* ATAN.C
  339.  * Calculates arctan(x)
  340.  * Range: -infinite <= x <= infinite (Output -pi/2 to +pi/2)
  341.  * Precision: +/- .000,000,04
  342.  * Header: math.h
  343.  * Author: Max R. Drsteler 9/15/83
  344.  */
  345.  
  346. double atan(x)
  347. double x;
  348. {
  349.     double xi, q, q2, y;
  350.     int sign;
  351.  
  352.     xi = (x < 0. ? -x : x);
  353.     q = (xi - 1.0) / (xi + 1.0); q2 = q * q;
  354.     y = ((((((( - .0040540580 * q2
  355.              + .0218612286) * q2
  356.              - .0559098861) * q2
  357.              + .0964200441) * q2
  358.              - .1390853351) * q2
  359.              + .1994653599) * q2
  360.              - .3332985605) * q2
  361.              + .9999993329) * q + 0.785398163397;
  362.     return(x < 0. ? -y: y);
  363. }
  364.  
  365.  
  366. /* FLOOR.C
  367.  * Returns largest integer not greater than x
  368.  * Author: Max R. Drsteler 9/26/83
  369.  */
  370. double floor(x)
  371. double x;
  372. {
  373.     double y;
  374.     int ix;
  375.  
  376.     y = 0.0;
  377.     while (x >= 32768.0) {
  378.         y += 32768.0;
  379.         x -= 32768.0;
  380.     }
  381.     while (x <= -32768.0) {
  382.         y -= 32768.0;
  383.         x += 32768.0;
  384.     }
  385.     if (x > 0.0) ix = (int) x;
  386.     else ix = (int)(x - 0.9999999999);
  387.     return( y + (double) ix);
  388. }
  389.  
  390.  
  391. /* CEIL.C
  392.  * Returns smallest integer not less than x
  393.  * Author: Max R. Drsteler 9/26/83
  394.  */
  395. double ceil(x)
  396. double x;
  397. {
  398.     double y;
  399.     int ix;
  400.  
  401.     y = 0.0;
  402.     while (x >= 32768.0) {
  403.         y += 32768.0;
  404.         x -= 32768.0;
  405.     }
  406.     while (x <= -32768.0) {
  407.         y -= 32768.0;
  408.         x += 32768.0;
  409.     }
  410.     if (x > 0.0) ix = (int) (x + 0.999999999999999);
  411.     else ix = (int) x;
  412.     return( y + (double) ix);
  413. }
  414.  
  415. /* EXP.C
  416.  * Calculates exponens of x to base e
  417.  * Range: +/- exp(88)
  418.  * Precision: +/- .000,000,000,1
  419.  * Author: Max R. Drsteler, 9/20/83
  420.  */
  421.  double exp(xi)
  422.  double xi;
  423.  {
  424.     double y, ex, px, nn, ds, in;
  425.  
  426.     if (xi > 88.0) return(1.7014117331926443e38);
  427.     if (xi < -88.0) return(0.0);
  428.     ex = 1.0;
  429.     while( xi > 1.0) {
  430.         ex *= 2.718281828459; /* const. e */
  431.         xi--;
  432.     }
  433.     while( xi < -1.0) {
  434.         ex *= .367879441171;  /* 1/e */
  435.         xi++;
  436.     }
  437. /* Slow, but more precise Taylor expansion series */
  438.     y = ds = 1.0; nn = 0.0;
  439.     while ((ds < 0.0 ? -ds : ds) > 0.000000000001) {
  440.         px = xi/++nn;          /* above precision required */
  441.         ds *= px;
  442.         y += ds;
  443.     }
  444.        y *= ex;
  445.  
  446. /*  Chebyshev polynomials: fast, but less precise then expected!
  447.  *    xi = -xi;
  448.  *    y = (((((.0000006906  * xi
  449.  *        +.0000054302) * xi
  450.  *        +.0001715620) * xi
  451.  *        +.0025913712) * xi
  452.  *        +.0312575832) * xi
  453.  *        +.2499986842) * xi
  454.  *        +1.0;
  455.  *     y = ex / (y * y * y * y);
  456.  */
  457.     return(y);
  458. }
  459.  
  460.  
  461. /* LOG10.C
  462.  * Approximation for logarithm of basis 10
  463.  * Range 0 < x < 1e+38
  464.  * Precision: +/- 0.000,000,1
  465.  * Header: math.h
  466.  * Author: Max R. Drsteler 9/15/83
  467.  */
  468.  
  469. /* Method of Chebyshev polynomials */
  470. /* C. Hastings, jr. 1955 */
  471.  
  472. double log10(x)
  473. double x;
  474. {
  475.     double xi, y, q, q2;
  476.     int ex;
  477.  
  478.     if (x <= 0.0) return(0.); /* Error!! */
  479.     ex = 0.0; xi = x;
  480.     while (xi < 1.0 ) {
  481.         xi *= 10.0;
  482.         ex--;
  483.     }
  484.     while (xi > 10.0) {
  485.         xi *= 0.1;
  486.         ex++;
  487.     }
  488.     q = (xi - 3.16227766) / (xi + 3.16227766); q2 = q * q;
  489.     y = ((((.191337714 * q2
  490.           + .094376476) * q2
  491.           + .177522071) * q2
  492.           + .289335524) * q2
  493.           + .868591718) * q + .5;
  494.     y += (double) ex;
  495.     return(y);
  496. }
  497.  
  498. /* PW10.C
  499.  * Calculates 10 power x
  500.  * Range: 0 <= x <= 1
  501.  * Precision: +/- 0.000,000,005
  502.  * Header: math.lib
  503.  * Author: Max R. Drsteler 9/15/83
  504.  */
  505.  
  506. double pw10(x)
  507. double x;
  508. {
  509.     double xi,ex,y;
  510.  
  511.     if (x > 38.0) return(1.7014117331926443e38);
  512.     if (x < -38.0) return(0.0);
  513.     xi = x; ex = 1.0;
  514.     while (xi > 1.0) {
  515.         ex *= 10.0;
  516.         xi--;
  517.     }
  518.     while (xi < 0.0) {
  519.         ex *= 0.1;
  520.         xi++;
  521.     }
  522.     y = ((((((.00093264267    * xi
  523.         + .00255491796) * xi
  524.         + .01742111988) * xi
  525.         + .07295173666) * xi
  526.         + .25439357484) * xi
  527.         + .66273088429) * xi
  528.         +1.15129277603) * xi + 1.0;
  529.     y *= y;
  530.     return(y*ex);
  531. }
  532.  
  533.